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1 Introduction 



In recent years there has been much interest in the properties of low-dimen- 
sional quantum systems. It seems to be worthwhile, therefore, to present 
a new algorithm to efficiently and accurately calculate a very large number 
of time-dependent pair correlations in the bulk for the premier example of 
such systems, namely for the transverse Ising chain0 Many results have been 
derived before for this modeljf] but the foremost method is still to calculate 
very large determinants for each data point, using either the appropriate 
determinant for the open chain with both spins far from the boundary [10] , or 
the infinite McCoy-Barouch- Abraham determinant [IT] for the closed chain. 

In this paper we shall use a set of coupled nonlinear differential-difference 
equations derived by one of us [121 [131, related to similar identities for the 
planar Ising model [T^ [15] . From these we can obtain multiple millions of 
data points in one single run. In spite of what some colleagues have told us — 
this work is partly an answer to their unbelief — our nonlinear equations are 
highly effective for numerical computations as we shall demonstrate. We shall 
concentrate on the zero-temperature casejf] for which the initial conditions 
coincide with the diagonal correlations in the two-dimensional Ising model. 

Therefore, we shall first discuss in section 2 how to accurately calculate 
these diagonal correlations. For smaller separations we can use the algorithms 
of Jimbo and Miwa [12] or of Witte [13]. For larger separations we need 
several more terms in the old asymptotic expansion of T.T. Wu [1^ B5] . 
which we shall derive as a new result. 

In section 3 we introduce the transverse-field Ising chain. We shall also 
present new results for the asymptotic expansions for its correlations as a 
function of time. In section 4 we shall give details of how we solved the 
correlations numerically and show how well the numerical results agree with 

"'^The transverse Ising chain [T] is a special case of the XY model first introduced by 
Nambu [2] in 1950 for the isotropic zero-field case and generalized by Lieb, Schultz and 
Mattis [3] and Katsura [4J. 

^ We shall not consider the time-dependent zz correlations [3 [H [71 [H , which being 
only 2-by-2 determinants do not involve " Jordan- Wigner strings" of fermion operators. 

■^The case of finite temperature should work out fine also according to a calculation 
we have done with approximate initial conditions, utilizing the exponential decay with 
separation in the initial conditions. The case of infinite temperature has been studied 
before in great detail [Tni[ni[IHl[Il[2Dl[5Tl[121[13]. We shall use some of the older results 
for zero temperature [ni[ll[lSl[151[ll[lHl[2Sl[3ni[3I]- Finally, there are several partial 
results for finite temperature worth mentioning, e.g. [5^ [551 [Ml [551 [551 [571 [551 [5^1 liDlliT] . 
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the asymptotic expansions derived in section 3. We conclude with a few 
remarks in section 5. 



2 Results for the Two-Dimensional Ising 
Model 

In this section we shall present results for the diagonal pair correlations of 
the two-dimensional Ising model on the square lattice and its dual model. 
These models are described by the reduced interaction energies (with factor 
P = l/k^T absorbed) 

), (1) 

m n 

and 

1 '^m,n'^m+l,n (2) 

m n 

It is convenient to introduce the following short-hand notations for the in- 
teraction parameters: 

ti = tanh Ki = e'^^i , = tanh K2 = e'^^i , (3) 

and 

^1 = sinh(2fs:i) = l/sinh(2is:*) = 1/S;, 

S2 = sinh(2ir2) = 1/ sinh{2Kl) = l/Sl (4) 
We shall restrict the elliptic modulus 

k = S1S2 (5) 
to < /c < 1, as the case k > 1 is included via the dual model0 defining 

C{m, n) = (o-o,oO-m,n) = E c^o.oo-m.ne"^^ / ^ e"'^^, (6) 

{a} / {a} 

C*{m,n) = (cro,ot^m,n)* = E ^o,ot^m,ne"^^* / E e"'^^* ■ (7) 

{a} / {a} 



*The case A; < follows by simple gauge symmetry, but will not be used. 
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As already said before, for the initial conditions of the quantum chain, we 
need to calculate C{m,n) and C*{m,n) for m = n. This can be done using 
the quadratic recurrence relations of Jimbo and Miwa ^2] or Witte [13] 

When T 7^ Tc, we can use the well-known Toeplitz determinants [I5l H6] . 
It is, however, more efficient to use the equations provided by Jimbo and 
Miwa [12]. Denoting C*{n, n) = An and C{n, n) = C„, these can be rewritten 
as: 



= + k-'A-B-)/i{2n + 3) A J, (8) 

^^+l = (^n+l^^ + Cni?n)M„, (lO) 

= + C-^^D-^,)/ii2n + 1)AJ, (11) 

^n+l — i^n+l^n ~ ^n+l^n+l) / ^ (12) 

B^^, = (k^'A^^^.B^ + B^^,D^^,)/A^ , (13) 

which can be solved iteratively, in the above order for n = 0, 1, ■ ■ ■, starting 
from the initial conditions 

A,=B,=C,= 1, A, = 2E(A;)/7r, 

B+ = D, = k', C+ = = 1/A;', D+ = C, = 0, (15) 
A+ = 2(2 E{k) - K{k))/{nk'), Bq = 2k'{K{k) - E(A;)) /vr. 



where K{k) and E(/c) are the usual complete elliptic integrals, k' = Vl — k'^. 
From these equations one can calculate systematically more A^s and C„'s 
by iteration, keeping sufficiently many digits, until the asymptotic regime is 
reached. 

The required asymptotic expansions can be obtained from the Painleve 
VI equation (P VI) of Jimbo and Miwa [12] , extending the expansions of Wu 
[H] — or also eqs. (2.46) and (3.27) of Chapter XI of [15] — beyond second 
order. For solving PVI iteratively, we need information from the leading 
terms in the low- and high-temperature expansions of the connected pair 
correlations. 

The high-temperature expansion of C{m,n) is well-known and is often 
treated in graduate courses. Without loss of generality we can take m,n > 
and give the leading high-temperature term as 

c,„.„).(!!i±^o?«^^i::f. (16) 

ml nl ml nl 2™-+"- 
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where the combinatorial factor counts the number of staircase walks from 
(0, 0) to (m, n) on the square lattice. For the leading term in the low-tempera- 
ture expansion of the connected pair correlation 

C:{m,n) = C*{m,n)-{l-k'Y/\ k^S,S, = -^, (17) 

subtracting the square of the spontaneous magnetization, we can use (10a), 
(10b) or (10c) of [15]. When doing so we need to replace M ^ n, N ^ m, 
as currently it is more common to use the horizontal coordinate as the first 
one. Rewriting eq. (10b) of [T5] as 

(C*{m - 1, n) C*{m + l,n)- C*{m, nf) 

+S^{C{m,n - 1) C{m,n + 1) - C{m,nf) = 0, (18) 

we can ignore the contributions of C*(m, n) and C*{m + l,n) and use A; — > 0. 
We arrive at 

C*{m - 1, n) ^ S^(C{m, n)^ - C{m, n - 1) C{m, n + 1)) , (19) 
or equivalently 

C*{m, n) ^ S^{C{m + 1, n)^ - C{m + 1, n - 1) C{m + 1, n + 1)) , (20) 

with the leading order solution 

^ {m + n)l{m + n + iy. S^-^^Sj-^^ 
' ' ^ m\{m + l)\n\{n + l)\ 22™+2"+2 " ^ > 

Here the combinatorial coefficient can be recognized as a Narayana number, 
as it counts all staircase polygons needed in the leading order of the low- 
temperature expansionlf] Setting m = n, we get|^ 

X (2^)! k'' , (2n)! (2n + 1)! A;2"+2 ^ ^ 



^We thank Dr. Xavier Viennot for pointing this out at the Dunk Island Conference of 
2005. The g-analogues of the Narayana numbers were aheady studied by MacMahon four 
decades before Narayana rediscovered them, see [47] and references cited therein. 

^These results have been advocated by Dr. Ranjan Kumar Ghosh [15] as initial condi- 
tions for the derivation of the first few terms in the high- and low-temperature expansions 
using PVI, with the second result in presented as a conjecture. However, there is a 
much more efficient way to derive such expansions [491 150j . 
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Following Jimbo and Miwa, we define 



t = — , (T„ = t(t - 1)— lnC(n,n) - -t 



A;2 



dt 



< = t(t-l)^lnC*(n,n)-^, 
and then (T„ and cr* both satisfy the PVI equation 



(23) 



t{t 



d^cr„ 



dt2 



— n 



(t - 1)^: Cfr. 



+4- 



d(Tn 
dt 



dt 



dt 



^dOn 
' dt 



Substituting 



C(n, n) 



^-n/2 



m Vj/2\ j-2s 

exp ( E E 



j=l s=0 



/^(l - l/t)V4 

for increasing values of for m = 1, 2, . . ., while defining 

t + 1 1 + ^2 



(24) 



(25) 



X 



t 



(26) 



(x = Xg on page 256 of [35J]), we can easily solve all the Pj,s's. Whenever m 
is even, we also need fl22l) in the limit or a; 1 in order to extract 
Pm,m/2- Similarly, substituting 



-n-l 



C:{n,n) 



27™2(1 - 1/t)' 



exp ( E E 



(27) 



for increasing values of m = 1,2,..., with the same x {x = —x'^ in [15] now), 
we obtain the P*jjs using (12^ in addition for the whenever m is even. 

Using short Maple programs we went up to m = 20 and ended with 



_ 74074237647505 

^20,10 ~ oooocno ' ^'20,10 



673835095036826977 



(28) 



8388608 ' "^'^'''^ 10485760 

We have used all these results, but here we only give the results up to m = 10, 
i.e. 

A;" ( X {x^ - 1) x(25x2 - 27) 

' 8^ ^ 16n2 384n3 



C(n, n) 



^(1 - A;2)V4 



exp 



+ 



(x^-l)(13a:^-5) _ x(1073x^ - 1830x^ + 765) 
^ 128^4 5120n5 

(a;2-l)(412x4-425a;2 + 61) 
768n6 

x(375733x^ - 886725x^ + 660723x^ - 150003) 

229376n7 

- l)(23797x6 - milx"^ + 18055^2 - 1385) 

4096^8 

a;(55384775a;^ - 167281524a;6 + 179965314x4 

-79479684x2 + 11415087) 

2359296n9 

(x^ - l)(2180461x« - 5127404a;^ + 3945946^^ 

-1048244x2 + 50521) , 

+ 20480nio '"' ^ ^ 



+ 



and 



fc2"+2 / 7x 17x2- 10 901x=^- 783x 

= t; — TT, 7^ exp - — + 



+ 



27^2(1 - P)2 ^ V 4n 8n2 192n3 
899x4 _ 1062x2 + 194 _ 131411x5 - 196770x3 + 66375x 

64n4 2560n5 
8359 Ix^ - 151767x^ + 75033x2 _ g^gg 
^ 384n6 

17052139.r' - 36116187r^ + 23770797r'^ - I102125x 

16384n7 

11282939xS - 27723492x6 + 22515930x4 - 6419700x2 + 344834 

^ 2048^8 

37620804281x'^ - 104587369452x^ + 101707083486x5 

-39418182684x3 + 4677930225x 

1179648^9 

2049064082x^0 - 6360721245x^ + 7210080180x6 

-3544939170x4 + 670637250x2 - 24119050 \ /on^ 
+ 10240;^^^^ "')■ ^'^^> 

At this point we have an effective way of calculating the connected di- 
agonal correlation functions of the square-lattice Ising model and their log- 
arithms, by iteration up to a certain distance using many-digit precision as 
this procedure is unstable, and by asymptotic expansion for larger distances. 
This way we also have the zero-time zero-temperature xx-correlations of the 
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transverse Ising chain for < A; < 1. As we will need very large numbers 
of log C(n,n) and log C*(n,n) at fixed k we will have to do the sums over s 
only once saving computation time. 

For the self-dual case = 1, with C*{n,n) = C*{n,n) = C{n,n) there 
is a simple formula, which was already known to Onsager [51], discussed in 
detail by McCoy and Wu [IS], and implemented by us as: 

logC(0,0) = 0, logC(n,n) = logC(n-l,n-l) + r„, (31) 

n = log(2/7r), r„+i = r„ - log (1 - 1/(4^2)). (32) 

For this case we also have the large-n asymptotic expansion [SJ H^l SS] , but 
we did not need it for this paper. 

When k is close to 1, the results are described by Painleve III or V (PHI 
or P V) scaling limit results [l2l HHl |52] . Vaidya and Tracy [25j have analyzed 
the time-dependent pair correlations of the transverse Ising chain in this 
scaling limit. McCoy, Perk and Shrock [2Sl[2Z| found the Painleve II scaling 
limit describing the crossover from the space-like to the time-like regime for 
this quantum chain model. 



3 Quantum Ising Chain 

From now on we shall consider equilibrium bulk properties of the ferromag- 
netic Ising chain in a positive transverse field. The sites are labeled by 
integers j, (with —A/" < J < A/", and A/" ^ oo in the thermodynamic limit). 
We shall take the usual spin-i operator basis for the quantum chain 

j-th 

(33) 
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The spin operators obey the usual Schrodinger time-dependence in the 
Heisenberg picture 



af{t) = e'^Vj' e"^^* (in units for which k = 1, a = x,y, z), (34) 
with Hamiltonian 7i (or dual Hamiltonian Ti.*) 



oo 



^ = -I E iJ^]^]+i + Bcr]), H* = ~\ E (i?aja^:,, + Ja|), (35) 

j = ~-OD j = — 00 

with J, B > 0. The dual chain corresponds to the interchange of J and 
The pair correlation function 



Tr (e'*^aj e-"^aj^„ e"^^l 
Tr (e-f^'H) 



X.it) ^ {a]it)a^,J ^ (36) 



and the dual X*{t), with H replaced by H*, satisfy [El [13] 

' x„(t)x„(t) - x„(t)2 = 52(x:_i(t)x:^i(t) - x:itf) , 

< 

. x:{t)x:{t) - x:xtf = j2(x„_i(t)x„+i(t) - x„(t) 



(37) 



for all inverse temperatures (3 = 1/k-QT. At the critical field B = J this 
reduces to 



X„(t)X„(t) - X^tf = j\Xn-l{t)Xn+l{t) - Xnitf). (38) 

For real times t these equations represent a discrete nonlinear generaliza- 
tion of a system of hyperbolic partial differential equations. Therefore, we 
expect the initial-value problem to be stable and this is indeed the case for 
sufficiently small time steps. On the other hand, for imaginary times t, as 
needed for the calculation of the susceptibility, we have an elliptic system and 
we should treat these equations as a boundary value problem or work with 
an increasingly high number of digits as \t\ increases, as was done before 
for the (euclidian) two-dimensional Ising model, first for the uniform case 



^Strictly spoken, the dual positions are at the half-integers, but we can relabel them to 
the integers. Also, as the bulk thermodynamic limit for pair correlations has been well- 
established and is independent of boundary conditions, we have formally replaced A/" by oo 
in the bounds of the summations in (j35|) . This thermodynamic limit is to be understood 
to be part of the pair correlation's definition in the following. 
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[531 IMl ESI ES] and later for quasiperiodic cases In fact, the 

susceptibility for the transverse Ising model has been calculated in the mid 
1980s as an extremely anisotropic limit of the corresponding 2d Ising result 
[531 [Ml EH E6]. 

Before describing the algorithm used, we shall first derive some long-time 
asymptotic expansions, extending the results on the last page of [TT], while 
also implicitly correcting several misprints there. We start with adopting the 
results (2.1) and (2.7) in [25] to our case For T = 0, 



XRit) 



exp 



n=l 



where 



1-^) X(--)exp 



_J2p{2n)y (0<5<J), (39) 

oo 

(0< J<S), (40) 



n=l 



1 

2n 



211 



i=i 



(41) 



and 



x(2«-i)(i?,t) = r d0i 

2n-2 



d0: 



2n-l 



i=i 



L(02„_i)M(</.2„-i,0i). (42) 



Here, 



exp — iR(j) — itA(0) 



M(0,0') = 



47rA(0) 
A(0) - A(0') 



AJB sin 



sin 

(1(0 + 00) A(0) + A(0O 



(43) 
(44) 



^These latter cases could be done because the general quadratic identities of [12] apply 
also to Z-invariant cases [61] with further implications for correlations and susceptibilities 

[MillMlISl]. 

^There are misprints in (1.2) and (1.3) of _25J, which differ from the actual definitions 
used in [TTJ [T51 [501 HH ISSj ■ The difference means that the results change to their complex 
conjugates, or equivalently to the replacement t —t. 
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(45) 



X{(j)) = J + -2 JB cos (f). (46) 



It is easily checked that the integrands in fl4Tl) and fl42|) are periodic modulo 
2tt in all the (pj. Therefore, before applying the stationary phase method for 
large t, we can shift the integration bounds and consider the stationary phase 
points (pj = and (pj = ir = —ir as internal point for the integration over (pj. 
Expanding the integrands of ( HTi) and ( H2i) in power series both in (pj and in 
(pj — 7T for all j , keeping only the quadratic terms in the exponentials of the 
L(0j)'s, we can then use 

d0e±^'^'^V'""' = 0, d0e±'"<^V = e±(2"+i)W4 ^i^+p _ (47) 

J —00 J— 00 a"^2 

for a > 0. The leading terms come from X^^^ and F*^^^ and to leading order 
one finds 

as the stationary phase points contributing are alternatingly and vr, starting 
with 01 = or vr. In higher orders of 1/t more and more non-alternating 
combinations come in and it soon becomes very tedious and unpresentable. 
We shall, therefore, only give the results for i? = in (l39i) and (110|) here. 
We use the abbreviations 

B J 

k = -, t = Jt (for B <J), k = -, t = Bt (for B > J). (49) 
J B 

Then, for B < J and t oo, 

1/4 / /ce~^** 



Xoit) = il-k 



,2 



+ 



V 27r Vl - k'^t 
ik{A - 3e)e-^'' e-2i(i-fc)i ^-2^(1+^)? 



+ 



+ 



07r(l-P)^/'t' 87r(l-A;)'t' 87r(l + A;)'t' 
(4 - 6OP + 80A;^ - 33fc6)e-2^* 

647rA;(l -P)^/'t^ 
i(l -9k + k'^)e-m-k)t + 9A; + e)e-n^+k)t 



32nk{l-kft^ 32nk{l + kft^ 
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327r2(l-P)'t^ 



(50) 



whereas, for B > J and t —>■ oo, 



A'o(«) 



-i(l-k)t-m/4: 



^-i(l+fe)t+7ri/4 



y V2^(A;(1 - k)t ) 
(1 -3A; + fc2)e-*(i-^)*+"*/^ 
8\/2^(fc(l - k)t 



1/2 / x-\l/2 



27r + A;)t 



3/2 



8\/2^(fc(l + k)t 



3/2 



3 (3 - lOA; + 17P - iQ/c^ + 3fc^) e^ 

- \ 5/2 



-•t(l-fc)t-7rj/4 



128^2^ (A;(1 - k)t 
3 (3 + lOA; + 17A;2 + lOA;^ + 2>k'^) Q-Ki+k)t+^i/i 



- \5/2 



128v^(A;(1 + k)t 

^3/2 g-i(3-fc)f+7rj/4 

4(27r)^/'(l- A;)'(l + A;)'/' t^'^ 

^3/2 g-i(3+fc)f-7ri/4 



4(27r)^/'(l + A;)'(l-A;)'^" t 



1/2 r 5/2 



+ ■ 



(51) 



The lowest order terms that decay as and have been given before 
in eq. (5.14) of [TT], and they agree after one corrects some misprints in [TT]. 

At the critical field B = J, we can use the Painleve V equation of [26] to 
obtain many terms in the long-time asymptotic expansion of Xo(t). Writing 



Xo{t) = e'''/Vo(x), aoix) = x 
we must solve 



dlogro(x) 
dx 



X = 2iJt, 



'x—Y A(x—-a 
, dx'^ I \ dx 



da 
dx 



X- 



'da' 
dx , 



a 



0, 



(52) 



(53) 



requiring the large-x asymptotic behaviors 



^2 N Li"+lJ ^n+2-2m 

^ n=0 m=0 



(54) 
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^2 A ^ '"^"^ ^n-2m 

logro(x) - -- + log^^y^ + E E ^"^^^ (55) 



and 



where ^ ^ 

(Here we use the more common real definition of A without the e"""^/^ factor 
of [26]. This factor is here supphed by the i^l^ of a:^/^.) 

We can solve the coefficients Cn,m recursively from the sets of linear equa- 
tions obtained by using the leading term of order 0(x'-^~^''/^) when substi- 
tuting (with = 0, 1, . . ., in succession) into (1551) . After each step in 
this process there are two undetermined coefficients left that are determined 
in the next two steps. The coefficients ,„ follow recursively for n = 3, 4, . . ., 
from ^ 

(n - 2m) dn,m = -Cn-2,m - -{u - 2) dn-2,m-l, (57) 

using also 

dn,0 = ) (^2 1 = 0, d2n,n = —- — ■ (58) 

n n 
The final result, after exponentiating logro(x), is 

, , A z 9z 2z'^-l 297 z 



{x/2y/^\ xV2 8a;3/2 8x^ 2'^ x^/^ 

15 z^ 7587 z 489^2- 81 1024 _ ^027035 ^ 

9387 ^2 76800 - 43594695 z 

+ 



29 2i8a;iV2 
851427 z^ - 90072 9094144 z^ - 4418168445 z 



+ 



+ 



213 2,6 222a;13/2 

22520925 529007616 - 260700970635 z 



215 2,7 225x1^2 

768 - 1368815805 z^ + 108135000 



218 
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+ 



+ 



258931316736 23 - 139999291654995 ^ 

231 3,17/2 

59904^^-47100085335^2 
220x9 

17114655467520 - 10543684346529075 z 

234 a;19/2 

24846336 - 14491877193315 + 908002224000 
2469862452602880^3 - 1758124895330287575 ^ 

238 2;21/2 

1113772032 2^ - 616117763829645^2 
_ 195069658230835200^3 _ 160841975585736493125 ^ 

241 2;23/2 

+ 0(x-i2)j (59) 

The first line of tlie above equation agrees witli (43) in [26]. Here we liave 
given a more efficient method to extend the asymptotic expansion, using real 
rational coefficients onlyS 

4 Numerical Results for Quantum Ising 
Chain 

At zero temperature, the initial values relate to the diagonal correlation 
function (crooO"„„) of the two-dimensional Ising model discussed in detail in 
section 2. More precisely, 

X„(0) = (aooa„„), X„(0) = Xo(0) 5„o = -iB {a^) 5no- (60) 



Here, k = J/ B = S1S2 < 1 corresponds to the 2d Ising model with T > 
while the case k = B / J = SIS2 < 1 corresponds to T < Tc. The initial value 



^''Replacing the term z^/x in (|54[) by —4Kx + (z + Kj z)\fx and summing m from to 
n + 2, we can easily extend also (39) in [26j . correcting two misprints in the highest order 
terms there in the process. One must then identify K — ~iah and z — be~^'^~^^^'^ . 
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of the first time-derivative X„(0) vanishes for ri 7^ and for n = it is given 
by the magnetization of the transverse Ising chain [H HI 

{al) = -m. ^=^^1' (61) 

vr B 

K) = ^(E(fc)-(l-A:')K(fc)), fc = f<l, (62) 

where K(A;) and E(fc) are the complete elhptic integrals of the first and second 
kind. For J = B, we have the simple results, see e.g. [26] , 

^"(0) = [-) n (1 - 4^') ' ^n(O) = - ^nO. (63) 

We have tacitly assumed n > 1 in the above. But, because of reflection 
symmetry and reality of the Hamiltonian Ti, we have 

X_„(t)=X„(t)=X„(-t)*, (64) 

where the asterisk denotes complex conjugation. 

The results in (!63l) are sufficient by themselves for B = J. But, also for 
B ^ J we know the initial conditions now to high precision from earlier work 
and the newly derived asymptotic expansions for large n in section 2. From 
these results we can find as many time-derivatives at t = as we want using 
the differential equations fl37j) and fl38l) . We can then calculate X„(5t) and 
Xn{St) and their duals if -B 7^ J, to high precision using Taylor expansion to 
five orders for sufficiently small 6t, or even more orders. Repeating this pro- 
cess N times we can calculate Xn{N6t) and Xn{N6t) (and their duals) using 
initial conditions in the "past light-cone." This is a discrete version of the 
method of characteristics. This is appropriate, since (1371) and (l38l) are non- 
linear differential-difference equations of hyperbolic type. In the continuum 
limit and for B = J, logX„(t) satisfies the two-dimensional Klein-Gordon 
equation [13]. In the more general continuum scaling limit B J we obtain 
the hyperbohc sine-Gordon equation [T31 125] . 

4.1 The case of critical transverse field B = J 

When B = J, me must solve fl38|) with initial conditions fl63|) . It is easily 
checked that X„(t) only depends on n and t = Jt. Equivalently we can 
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choose the unit of time such that J = 1 and t = t, which we shall do in this 
subsection. 

It is also convenient to work with the new variables 

Ut)^logX„(t), Ut)^^ = §^y (65) 
r]nit) = "^J. = exp(^„+i + - 2^n)- (66) 



Knowing these one can get the higher time-derivatives from fl38|) . i.e., 

Ut) = Vn-l, (67) 

Ut) = Vn{in+l+in^l-2in), (68) 

Ut) = + 24 + (4+1 + 4-1-24)'), (69) 

'in{t) = Tjni^ 4+1 ~^ 4-1 ~2 4 

+3(Cn+l + in-l — 2^n)(^n+l + ^n-1 — 2^„) 

+(4+1 + 4-1-24)'), (70) 



and we can then calculate ^„(t + 5t) and ^nii + ^i) in Taylor series in 5t to 
fifth, respectively fourth order. 

We have implemented this on a Macintosh computer using GNU Fortran 
with 5t = 10~^. Starting at t = we did = 300, 000 time steps, going up 
to t = 30. Therefore, we calculated X„(t) in the triangle in the n-t plane 
with the three corners {n,t) given by (±6x10^,0) and (0,30). Some of the 
results are given in Fig. 1. Comparing the Fortran output for Xo{t) with the 
asymptotic expansion ( l59i) . we concluded that we have already about 12 place 
agreement for t > 20, even though the Fortran program only worked with 
15 place floating point complex numbers. This is as good as we could have 
hoped for and it demonstrates the stability of the procedure. It is a strong 
indication that all the Fortran output within the "characteristic triangle" is 
equally accurate. 

Next, we can use the general identities of Lajzerowicz and Pfeuty [57] . 



1 dXnjt) 

B dt 

1 d'Xnit) 



Cr.{t) ^ (^J(i)^|+n) = -(^J(t)^J+„) = ^, (71) 



F„(t) = (-|(t)^J+n) = (72) 
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Figure 1: Xo(t) for i? = J = 1. On the left are plotted |Xo(t)| and argXo(t), 
while on the right ^Xo{t), and 53Xo(t) are shown. 

which are valid for all J and B and also for finite temperatures. These 
quantities are also calculated in the Fortran program with 5 = J = 1 as 

Cn{t) = Ut) Ynit) = -{Lit) + Utf) e«"W. (73) 

Their asymptotic expansions for n = follow immediately from flMJl) and 
agree with the Fortran data with a similar accuracy as the Xo{t) did. We 
have plotted Yo{t) in Fig. 2. 

Finally, we remark that we could have made the plots using only the 
Painleve V equation of [SHI [2Sl EH]- We would, however, not have been 
able to calculate the about 10^^ other very accurate data points within the 
"characteristic triangle" that way. 

4.2 The case of noncritical transverse field B ^ J 

In this subsection we assume that J < B in our transverse-field Ising chain, 
so that B* < J* for the dual Ising chain with B* = J and J* = B. We 
choose time units such that i? = J* = 10 whence J = B* = k and (l37|l is 

^^Equivalcntly, we may replace t by t = Bt = J*t, k = J/B in the following formulae. 
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Figure 2: lo(t) is plotted for B = J = 1. On the left ^Yo(t) is shown and on 
the right QXait). 

to be replaced by 

x^it)x^it) - x^itf = x:_,it)x:_,,it) - x:itf , 

(74) 

x:{t)x:{t) - x:{tf = fc2(x„_i(t)x„+i(t)-x„(t)2). 

At the critical point T = 0, B = J, we found it better to work with the 
logarithm ^„(t) of the pair correlation function Once we leave the 

critical point it becomes even essential to do so because of the exponential 
decay of the connected pair correlation with separation n, as we need initial 
conditions for very large n. Therefore, we introduce ^„(t) and ^*{t) via 

XrXt) = (l-A;2)V4e«.W, (75) 
X:{t) = (l-fc2)i/4(^^^cW), (76) 

Substituting these into flTil) these equations become 

^„(if:) = e^"-i"^^" +e^"+i~^^" -2e^""^^" 
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(78) 



These expressions were written such that none of the real parts of their 
exponents goes rapidly to +00. 

The initial conditions for ^n(t) and ^^(t) are determined by ( l29l) . ( l30l) . 
(!60|) . (!6T|) . and (l62l) . One immediately sees that ^n(O) and ^j^(O) roughly grow 
linearly with |n|, which is much better for the numerics than the exponential 
behavior of X„(0) and X*(0). Setting, as explained above, S = 1 in (160|1 for 
X„(0) and B* = k for X*(0), we find 



C(0) 



2iE(A;) 



TT 



2i [E{k) - (1 - P)K(A;) 



TT 1 



A;2) 



1/4 



5, 



nO- 



(79) 
(80) 



The higher derivatives at t = can then be evaluated using the pair of 
coupled differential equations ( 1771) and ( 1781) . 





Figure 3: Xo{t) for B = 0.7, J = 1: ^Xo{t) and '^Xo{t). 

We can now adopt the same strategy as for the critical-field case B = J. 
Given ^„(t) and ^j^(t) and their first five time-derivatives at a given t, we can 
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Figure 4: Xo(it) for B = 0.7, J = 1: ^Xo{t) and ^Xo{t). 

calculate + 6t) and + 6t) in a Taylor series through fifth order in 6t 
and C,n{t + St) and + in a Taylor series through fourth order. We have 
implemented this in a Maple program with 6t = 10~^ going up to t = 8.77. 

Our results for Xo(t) at = 0.7 are represented in Fig. 3, using both 
our Maple numerical integration results for t < 8.77 and our asymptotic 
expansions (I5T!) for t > 5.77. There is clearly excellent agreement between 
the numerical integration and the asymptotic expansiono 

We have also plotted Xoit) and lo(^) ioi t < 50 using the numerical 
results for t < 8.77 and the asymptotic expansions for t > 8.77. We have 
again used fl72|) to calculate Yo{t). Note that the lowest frequency, clearly 
seen in Fig. 4 for Xo(t) is hardly visible after two time derivatives in Fig. 5. 
The ratio of the two leading frequencies is 3/17, as follows from flSTl) . leading 
to a suppression factor of 9/289. 

Because of the very nature of the differential equations of \i2i used, we 
also got the results for the dual model from the same Maple run. They 
correspond to X„(t) at B/J = 1/k, or J/B = 0/7. We have plotted Xo(i) 
for t < 12 in Fig. 6. Again it is seen that the numerical results approach 
the asymptotic behavior well. Note that the vertical scales are smaller than 

^^We note that there is no agreement with the corresponding lowest-order asymptotic 
results in yjj because of several misprints there. 
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Figure 5: Yo{t) for B = .7 , J = 1: ?ftYo{t) and ^Yo{t). 

in Fig. 3, making the small differences more pronounced. We have also 
plotted Xn{t) and Yn{t) for t < 50 for the same values of the parameters. 
Unlike Fig. 4, there is now no pronounced two-frequency behavior, as other 
frequencies are down two orders in in (ISUl) . In addition, it is clearly seen 
that Xo(t) decays to the square of the order parameter (1 — fc^)^/^ 0.845, 
as it should [66] . 

5 Final Remarks 



The original motivation for this work came from two sides. First, we were 
told by some colleagues that solving the equations of [12] numerically leads 
to useless results inconsistent with earlier results, especially [H] . We believe 
that this point of view has been put to rest by the results in this paper. 
Secondly, during a brief visit in 2005 at the University of New South Wales 
one of us got to discuss an early version of |68|, [69] with his colleagues. It led 
us to calculate the small and large-i? expansions for the structure function 
(Sj '^0 '^j ) of Ising chain in a transverse field at zero temperature. This 
turned out to be nothing but the high- and low-temperature series of the 
"diagonal susceptibility" of the two-dimensional Ising model PJ2n{(^oo<^n,n), 



21 




Figure 7: Xo(t) for J = 0.7, 5 = 1: 5RXo(t) and QXo{t). 
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Figure 8: Yo{t) for J = .7 , B = 1: ^Yo{t) and ^Yo{t). 



the calculation of which is a small part of the calculation of ^\ ISO]. A 
different approach to the diagonal susceptibility in terms of form factors 
and ra-particle contributions has been given in [701 EI] revealing a lot of its 
mathematical structure. 

With the results of this paper it is now also very easy to calculate ex- 
tremely accurate values of the g-dependence of the structure function [681 ES] 
and of the second moment correlation length (as used e.g. in [72j). 

Furthermore, it is well-known that the correlation functions of the zero- 
field XY-model factor into two correlation functions of the Ising model in 
transverse field [201 122], just like correlations in the two-dimensional square- 
lattice dimer and fully-frustrated Ising models factor [731 Hi] ■ Therefore, our 
results directly apply also to those cases. 

Finally, as more and more results on time-correlations in the more general 
XXZ model have become available [7Sl[7Sl[77l[7Sl[7ni[Hni[SIl[S2],one may 
express the hope that some new powerful identities will be discovered also 
for this case that will facilitate numerical calculations dramatically. 
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